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Abstract 

Abelian sandpile models, both deterministic, such as the Bak, Tang, Wiesenfeld (BTW) model 
[P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987)], and stochastic, such as the 
Manna model [S.S. Manna, J. Phys. A 24, L363 (1991)], are studied on finite square lattices with 
open boundaries. The avalanche size distribution Pl(ji) is calculated for a range of system sizes, 
L. The first few moments of this distribution are evaluated numerically and their dependence on 
the system size is examined. The sandpile models are conservative in the sense that grains are 
conserved in the bulk and can leave the system only through the boundaries. It is shown that 
the conservation law provides an interesting connection between sandpile models and random walk 
models. Using this connection, it is shown that the average avalanche sizes, (n)i, for the BTW 
and the Manna models are equal to each other, and both are equal to the average path-length of 
a random walker starting from a random initial site on the same lattice of size L. This is in spite 
of the fact that sandpile models with deterministic (BTW) and stochastic (Manna) toppling rules 
exhibit different critical exponents, indicating that they belong to different universality classes. 
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I. INTRODUCTION 



Sandpile models have been studied extensively in the past fifteen years as a paradigm of 
self organized criticality (SOC) jl], 0, SOC provides a useful framework for the analysis 
of driven nonequilibrium systems which dynamically evolve into a critical state. At the 
critical state these systems exhibit avalanche dynamics with long-range spatial and temporal 
correlations, which resembles the behavior at equilibrium critical points. In sandpile models, 
defined on a lattice, grains are deposited randomly until the height at some site exceeds the 
threshold, thus becoming unstable. The unstable site is toppled and grains are distributed 
between its nearest neighbors, which may become unstable too, resulting in an avalanche. 
These models were found to be self-driven into a critical state in which the avalanche sizes 
follow a power-law distribution. The critical state, which can be characterized by various 
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171 Il8l Il9l |20j . These studies stimulated 
an effort to examine the utility of the SOC framework to the understanding of empirical 
phenomena such as earthquakes avalanches in granular flow and mass extinctions |2l| . 
To examine the dependence of the critical state on various properties of the models, dif- 



ferent sandpile models have been introduced. These include the stochastic model introduced 

and 



by Manna |22j. The issue of universality has been studied. Analytical studies \2c 
numerical simulations [25] indicated that the Manna model, which is stochastic, belongs to 
the universality class of the original model introduced by Bak, Tang and Wiesenfeld (BTW) 
which is deterministic (namely, has a deterministic toppling rule). However, numerical sim- 
ulations using an extended set of critical exponents showed that deterministic and stochastic 
models exhibit different scaling properties and thus belong to different universality classes 
26 . 27 . 28]. Further support for this result was obtained using multifractal analysis 
moment analysis [3^ as well as studies of sandpile models as closed systems j^l, 32, 33^ The 
crossover between the two classes was also studied . In the case of directed models it was 



shown analytically that deterministic and stochastic models belong to different universality 
classes (4, 



In this paper we present a connection between abelian sandpile models and random 
walkers on finite lattices, which is a result of the conservation laws. In the sandpile models 
each avalanche starts with the addition of one grain. The models are conservative in the 
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sense that grains are conserved in the bulk of the system and can leave it only through 
the boundaries. It is thus clear that under steady state conditions, the average number of 
grains leaving the system per avalanche is also one. Here, the avalanche size is defined as 
the number of hops of grains that are toppled from unstable sites during an avalanche. The 
avalanche size is thus equal to the number of toppling events (or relaxations) of unstable sites 
during the avalanche times the number of grains that topple in each event of this type. We 
show that in both the BTW and the Manna models, each grain moves like a random walker, 
starting at the site in which it was deposited, until it falls off the boundaries. Using these 
features we show that the average avalanche size is the same for the two models. Moreover, 
it is equal to the average path length of a random walker starting from a random site on 
the same lattice, until it falls off the edge. 

In order to demonstrate these properties we examine the avalanche size distribution of 
the Abelian sandpile models and the distribution of path length of a random walker on 
finite square lattices. The average path length of a random walker on a lattice of size L 
is calculated exactly using a method proposed by Walsh and Kozak The entire 

distribution of the path lengths of random walkers starting at random sites on the finite 
lattice is also calculated using a related method proposed by Soler [3^. The results are also 
compared to direct numerical simulation of the random walk. The avalanche size distribution 
of the Abelian sandpile models is obtained from direct numerical simulations as well as from 
an exact formula introduced by Dhar 

The paper is organized as follows. The sandpile models are described in Sec. II. The 
distributions of path lengths of random walkers on finite lattices and their averages are 
studied in Sec. III. The connection between sandpile models and random walkers is examined 
in Sec. IV. The simulations and results are given in Sec. V, followed by a discussion in Sec. 
VI and a summery in Sec. VII. 



II. SANDPILE MODELS 



Consider a sandpile model on a rf-dimensional cubic lattice of linear size L. Each site i 
is assigned a dynamic variable E(i) which represents some physical quantity such as energy, 
grain density, stress, etc. A configuration {E(i)} is called stable if for all sites E(i) < E c , 
where E c is a threshold value. The evolution between stable configurations is by the following 
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rules: (i) Adding energy. Given a stable configuration {E(j)} we select a site i at random 
and increase E(i) by some amount 5E. When an unstable configuration is reached rule (ii) 
is applied, (ii) Relaxation (or toppling) rule. If E(i) > E c , relaxation takes place and energy 
is distributed in the following way: 

E (i) -> £(i)-]TA£(e) 

e 

E(i + e) -> £(i + e) + A£(e), (I) 

where e are a set of vectors from the site i to some neighbors. As a result of the relaxation, 
E(i + e) for one or more of the neighbors may reach or exceed the threshold E c . The 
relaxation rule is then applied until a stable configuration is obtained. The resulting sequence 
of topplings is an avalanche which propagates through the lattice. 

Avalanches can be characterized by their size. The size s of an avalanche is the total 
number of toppling events that occurred during the course of the avalanche. In the models 
studied here the number of grains that topple from an unstable site is E c = J2 e AE(e). 
Throughout the rest of the paper we will denote the avalanche size by 

n = E c ■ s, (2) 

namely, by the number of hops of individual grains that take place during the avalanche. 
This will allow us to consider models with different values of E c on a common footing. 

The avalanche size distribution is denoted by Pl(ti), n = 0, 1, . . ., namely, the probability 
of a randomly chosen avalanche to be of size n. The normalization condition is given by 

oo 

$>L(n) = l. (3) 

n=0 

Numerical simulations show that the avalanche size distribution for a lattice of size L has 
the power-law form: 

P L (n)~n- TL , n = l,2,.... (4) 

For some of the models the results for tl exhibit a significant dependence on the system 
size. The critical exponents tl turn out to depend on the vector AE to be termed relaxation 
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vector. For a square lattice with relaxation to nearest neighbors it is of the form AE = 
(En, Ee, Es, Ew), where En, Ee, Eg and Ew are the amounts transferred to the northern, 
eastern, southern and western nearest neighbors respectively. The average avalanche size on 
a lattice of size L is given by 



The sandpile models considered in this paper are conservative in the sense that the grains 
are conserved in the bulk and can leave the system only through the open boundaries. When 
an avalanche reaches a boundary site, some energy is transferred out of the system (namely, 
dissipation takes place at the boundaries). The critical state is reached spontaneously in 
the limit in which the random addition of energy (or drive) is infinitely slow (practically it 
means that the next energy unit is added only after the previous avalanche is completed). 
This state is characterized by a power-law distribution of avalanche sizes [Eq. (J3J)]. In the 
critical state the added energy SE per avalanche, is balanced, on average, by the energy 
that flows out through the boundaries. Therefore, the average amount of energy leaving the 
system per avalanche is 5E. 

In the BTW model, E c = 4, SE = 1 and AE = (1,1,1,1). Since AE is a constant, 
this model is clearly deterministic. Note that since AE is independent of E(i), if an active 
site with E(i) > E c is toppled, it remains non-empty after the toppling event had occurred. 
A useful way to analyze the BTW model is by its toppling matrix A, which for an L x L 
lattice is a matrix of size L 2 x L 2 . Consider a pair of sites i = (i x , i y ), j = (j x ,j y ) and denote 
% = Li x + i y , j = Lj x + j y where < i x ,i y , j x , j y < L — 1. The matrix element Ajj where 
i 7^ j is the number of grains given to site j when site i topples (up to a minus sign). The 
number of grains leaving site i in such event is given by the diagonal element A^j. Therefore, 
the toppling matrix is : 



oo 



( u )l = nP L{n). 



(5) 



n=0 
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< — 1 i and j are nearest neighbor sites 



(6) 







otherwise. 



Consider the toppling of a given site i. It can be described by: 
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E(j) - E(i) - Ay-, (7) 

for all sites j. 

In the class of stochastic sandpile models, introduced by Manna, a set of neighbors is 
randomly chosen for relaxation 22| once a site becomes unstable. Such models can be 
specified by a set of relaxation vectors, each vector is assigned with a probability for its 
appearance. There are several models in this class. One of them is a two-state model with 
E c = 2 and two relaxation vectors (1,0,1,0) and (0,1,0,1) each one applied with a probabil- 



ity of 1/2 |40j. Another two-state model includes six relaxation vectors, namely (1,1,0,0), 
(1,0,0,1), (0,1,1,0), (0,0,1,1), (1,0,1,0) and (0,1,0,1), each one applied with a probability of 
1/6 [26] . In this paper we consider a two-state model in which each of the two grains of an 
unstable site is toppled randomly to one of the four neighbors (with probability 1/4 to each 
direction). There is no correlation between the directions picked for these two grains. The 
set of relaxation vectors includes all the 10 possible vectors AE of integer components for 
which E N + E E + E w + E s = 2. Each of the six vectors of the previous model appears with 
probability 2/16, while each of the four vectors: (2,0,0,0), (0,2,0,0), (0,0,2,0) and (0,0,0,2), 
appears with probability 1/16. This model will be called the unrestricted two-state Manna 
model. 

The average avalanche size (n) for the BTW model on a lattice of size L was calculated 
exactly by Dhar [6J. He showed that the matrix element A^j represents the average number 
of toppling events taking place at site j when a grain starts an avalanche after being deposited 
at site i. By summing all the elements of A -1 , using the eigenvectors and eigenvalues of A 
[dil l it was found that 

i x - cot 2 e k cot 2 gj 

W "L 2 (L + l) 2 ^sin 2 ^ + sin 2 ^' W 

where 

71771 (a\ 

m = 2(1+1)' [) 

for any integer m, and the summation over k, I is over all odd integers 1 < k, I < L. The 
dependence of (n) on the system size was found to be (n) ~ L 2 . This analysis was recently 



extended to dissipative abelian models |42j. 

III. RANDOM WALKS ON FINITE LATTICES 

Consider a random walker on a square lattice of size L. The walker starts at a random 
initial site i = (i x ,i y ), where < i x ,i y < L — 1. At each step the walker has four possible 
moves, to one of the sites = (i x +l, i y ), = (i x — l, iy), iu = (*x> S+l) an d *d = (ix, iy — !)> 
each picked with equal probability. The boundaries are open, and thus a random walker 
starting at any site (i x ,i y ) will eventually fall off the edge |37l. l38|. The number of moves it 
will make depends on the location of the initial site as well as on the particular realization of 
the random moves generating the path of the given walker. Therefore, there is a probability 
distribution Pi(n), where n = 1,2 ... , oo, that a random walker starting at site (i x , i y ) will 
fall off the edge after n moves. We will first calculate the average of this distribution given 
by 

oo 

( n i) = J2 n - Pi( n )- ( 10 ) 

71=1 

The boundary conditions are given by (ni) = for sites beyond the edge of the Lx L lattice, 
namely those for which i x = —1 or L, or i y = — 1 or L. Since on its way to the boundaries 
the walker must pass through at least one of the nearest neighbors of the site i, there is a 
relation between (ni) and the corresponding averages for its nearest neighbors of the form 

(ni) = \ «n iR ) + (n iL ) + (m D ) + (n, u » + 1. (11) 
This set of L 2 coupled linear equations can be written in a matrix form as 

^A(n) = 1, (12) 

where the matrix A is identical to the toppling matrix of the BTW model, given in Eq. (JBJ). 
The vector (n), consists of the L 2 components (n)j = (ni), where i = (i x , i y ) and % = Li x + i y 
(i = 0, 1, . . . , L 2 — 1). The L 2 -dimensional vector 1 is given by 1 = (1, 1, . . . , 1). 

In order to reduce the number of equations we shall use the symmetry properties of the 
square lattice, that has one horizontal, one vertical and two diagonal reflection axes. Any 
two sites with the same symmetry properties are called sites of the same "type" . Due to 
the symmetry it is sufficient to examine the sites in the triangle bounded by the vertical 
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axis from the center upwards and by the diagonal axis from the center to the upper - right 
corner (Fig. In this triangle there is one site of each type. The number of sites in the 
triangle is 



N 



(L + l)(L + 3) 



(13) 



when L is odd, and 



N 



L(L + 2) 



when L is even. The N linear equations for (n)i are of the form 



1 



D(n) = 1, 



where D is an N x N matrix. The matrix elements of D are: 



(14) 



(15) 



D 



'■j 



4-/(M) i = j 
-f(i,j) 



(16) 



where f(i,j) is the number of sites of type j that are nearest neighbors to a site of type i. 
For L = 3 Eq. ()15|) takes the form 



t 1 -1 \ / (n)i \ ( 1 



■3 1 



v -I 1 y 



(n>i 
(n) 2 

v («)3 y 



i 

v 1 / 



(17) 



and its solution can be easily found to be (n) = (4.5, 3.5, 2.75). Now, having explicit values 
for the (ri)j's, the average path length (n) of a random walker that starts at a random site 
on the 3x3 lattice is 



(n) 



(m) + 4(w 2 ) + 4(n 3 ) 
9 



3.277.. 



(18) 



The probability distribution pdn) of a walker starting at site i to fall off the edge after 
n moves can also be calculated [39j. One can then average Pi(n) over all lattice sites and 
obtain the probability distribution P/^(n), n = 1, 2, . . ., that a walker starting at a random 
site on the lattice will fall off the edge after n steps. This probability is given by 
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Pdn) = -Y.P^ n ). (19) 

i 

Note that for the random walk model, Pi(n) is defined only for n > 1, because the random 
walker must make at least one move in order to fall off the edge. The normalization condition 
will thus take the form 

oo 

£P z (n) = l. (20) 

n=l 

The moments of this distribution are given by 

oo 

(n*) L =J2n*.P L (n). (21) 

n=l 

The average path length (first moment) is given by Eq. (|21|) with q = 1, the second moment 
by q = 2 and so on. The calculation of Pi(n) is done recursively starting from the boundaries. 
The probability that a walker starting at site i will fall off the edge after n steps is given by 

Pi(n) = ^ bi R ( n _ 1 )+PiA n ~ 1 )+Pi v ( n ~ !) +Pi D ( n - !)] > ( 22 ) 

namely, it is the average over the four nearest neighbors, of the probabilities that a walker 
starting in one of them will fall off the edge after n — 1 steps. The boundary conditions are 
Pi( n ) — where i x = —1, L or i y = —1, L and n = 1, 2, . . . , oo, reflecting the fact that these 
indices represent sites that are already over the edge. The initial conditions for the recursive 
procedure for the calculation of Pi(n), n — 1, 2, . . ., are given by Pi(l) = 1/4 for all the edge 
sites, except for the corner sites, for which Pi(l) = 1/2. For all other sites Pi(l) = 0. The 
distribution Pi(n), n — 2, . . ., is calculated recursively using Eq. (|2^j) for all sites, starting at 
n = 2 and increasing n by 1 after each scan of the lattice. The average path length (n) of a 
random walker starting at a random site is then obtained from Eqs. (fT§|) and (|2*T|) . 

IV. THE CONNECTION BETWEEN SANDPILE MODELS AND RANDOM 
WALK MODELS 

Consider a grain deposited at a random site in the unrestricted two-state Manna model. 
It will initiate a small or a large avalanche and will typically stay on the lattice for many 
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subsequent avalanches until it will fall off the edge. Most of this time the grain is alone 
in its site. Whenever it will share the site with another grain, both of them will topple 
randomly (and independently) to nearest neighbor sites. The path of the grain on the 
attice is, in fact, a random walk since there is no correlation between one move to another 



43| . This path starts at the site into which the grain was deposited and ends at the edge 
site through which it leaves the system. The different walkers are uncorrelated since the 
directions chosen by the two grains that topple from an unstable site are independent. The 
only correlation is between the times that different walkers make their random moves. This 
temporal correlation appears because two walkers need to occupy the same site in order to 
move. 

This property, that the path of each grain in the sandpile model is a random walker, 
is not limited to the unrestricted two-state Manna model. It is a general property of the 
conservative-Abelian models, in which the grains are discrete entities. It is thus a common 
property of models that belong to both the deterministic and the stochastic universality 
classes. Consider, for example, the BTW model. It is convenient to consider the grains 
as distinguishable particles by naming each one of them according to the running number 
representing the order of their deposition into the system. When an unstable site topples, 
we will pick the grain that entered this site first, and choose randomly one of the four 
directions (N, E, S, W) for it to move. For the second grain we will choose randomly among 
the three remaining directions, and so on for the third and fourth grains. It is clear that 
there is no correlation between the directions of consecutive moves of each grain and no bias. 
Therefore, each grain follows a path of a random walker. Unlike the unrestricted two-state 
Manna model, in the BTW model there is correlation between the directions of different 
grains that topple from the same site, since they cannot move in the same direction. 

Each avalanche in the sandpile models starts with a new grain deposited randomly. There- 
fore, on average each avalanche drops one grain off the edge. Since the grains follow random 
walker paths from the random initial site to an edge site, the average number of hops that 
take place in a single avalanche must be equal to the average number of steps, that is re- 
quired for a random walker deposited randomly on the lattice to reach the edge. We thus 
conclude that the average avalanche size of a sandpile model on a lattice of size L is the 
same for the BTW and Manna models, and equal to the average path length of a random 
walker deposited randomly on the same lattice. 
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V. SIMULATIONS AND RESULTS 



To examine the connection between the random walk model and the sandpile models on 
finite lattices, we have performed numerical simulations of both systems. Direct simulations 
of the random walk model were performed on a square lattice of size L with open boundaries. 
In each run the walker started at a random site on the lattice. The random walk path was 
generated until the walker fell off the edge. The path length, namely the number, n, of moves 
it made from the initial site to the edge was recorded. From this data, the distribution of 
path lengths Pl{ti) was generated and its average (n) j, was calculated. The average (n)^ 
was also calculated using the Walsh-Kozak method 37|,|38|. The distribution Pt,(n) was also 
calculated using the Soler method j^J and (n) l was extracted from it. Direct simulations of 
the BTW and Manna models were performed, from which the avalanche size distributions 
were obtained. The average avalanche size for each model was calculated as a function 
of the lattice size L. The average avalanche size for the BTW model was also calculated 
using Dhar's formula [Eq. (JEJ)]. The values of (n) vs. L, obtained for both the random 
walk and sandpile models are shown in Fig. |21 They all coincide perfectly, except for some 
slight fluctuations in the direct simulation data for the larger values of L. This confirms 
the connection between the random walk and sandpile models. Fitting {ti)l vs. L to a 
polynomial function we obtain that (n) l = aL 2 + bL, where a = 0.14 and b = 0.56. 

While the averages are found to be the same for the sandpile models and the random 
walk models, the distributions Pl(tl) turn out to be different. To calculate the distribution 
Pi(rz) of path lengths of random walkers starting from random sites on a square lattice of 
size L, we used the Soler method 3^. The calculation was done for square lattices of sizes 
L = 32,64,128,256 and 512. To obtain a scaling function we rescaled n for each system 
size by the average path length (ti)l- The scaling function P{n/{n) L ) = (n) L ■ Pt(n) is 
shown in Fig. El on a double logarithmic scale. The scaling function exhibits a linear range, 
up to an upper cutoff around n/(n)L = 1- The slope of the linear range turns out to be 
approximately —1/2. This function can be considered in the framework of first passage 
problems of a random walks on a finite lattice Q|. 

The avalanche size distributions for the BTW and Manna models were obtained from 
direct numerical simulations for lattice sizes L = 32, 64, 128, 256 and 512. The rescaled 
distribution functions are shown in Fig. Eta) for the BTW model and in Fig. Efb) for the 
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Manna model. In both cases the data collapse is not complete, due to the finite size de- 
pendence of the critical exponent tl- Fitting the data for L = 512 to Eq. (0J) we obtain 



that r L = 1.12 ± 0.02 for the B' 
agreement with previous results 



W model and tl = 1-27 ± 0.02 for the Manna model, in 



22 



27 



The first three moments of the distribution of path-lengths of the random walk vs. L are 
shown in In Fig. |SJ The results were obtained both by direct simulation and by calculating 



Pi{n) using the Soler method 39j. Fitting these graphs to power laws in L we find that 
(n q ) L ~ L 2q (the slopes of the best linear fits are 1.98 ± 0.04, 3.97 ± 0.06, and 5.96 ± 0.08 
for q = 1, 2 and 3, respectively). 

For sandpile models, in the large system limit, the first moment of the avalanche size 
distribution scales like (n)i ~ L 2 . Higher moments are expected to scale like 

(n«) L ~L^\ (23) 

with a(q) > 2q for q > 2. The first three moments of the avalanches size distributions of 
the BTW and Manna models, vs. L, are shown in Fig. |H1 The slopes of the best linear fits 
for the BTW model are a(q) = 1.98 ± 0.02, 4.68 ± 0.04, and 7.52 ± 0.08 for q = 1, 2 and 3, 
respectively. The slopes of the best linear fits for the Manna model are a(q) = 1.97 ± 0.02, 
4.73 ± 0.04, and 7.48 ± 0.08 for q — 1, 2 and 3, respectively. These linear fits were obtained 
for lattice sizes in the range 64 < L < 1024. The results for the first moment are identical 
(within the error bars) for the two models, and coincide with the results for the random walk 
model, and thus confirm the conclusions of the analysis above. Surprisingly, the values of 
a(q) for the BTW and Manna models are approximately the same (within the error bars), 
also for q = 2 and 3. This is in spite of the fact that the avalanche size distributions of 
the two models are characterized by different exponents r^. This behavior has to do with 
deviations from the power-law behavior near the upper cutoffs of the distributions. The 
results for the higher moments are in agreement with those presented in Refs. jsnL l^. 

VI. DISCUSSION 

Power-law distributions were observed in a wide variety of natural systems as well as in 
economic systems, computer networks, linguistics and other fields. Some examples include 
the energy distribution between scales in turbulence the distribution of earthquake 
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magnitudes 1471 1 , th e distribution of city populations |48|, |49(, the distribution of income 



H 52,0 



HQ, 



the distribution of the number of links pointing to sites in the 



and wealth 

internet sj, 55| and the distribution of the frequency of appearance of words in texts |48(. A 
common feature of such system is that they consist of a large number of elementary degrees of 
freedom that interact with each other in a complex way. Power-law distributions typically 
appear when these interactions give rise to long-range correlations with no characteristic 
length-scale. 

Consider a power-law distribution of the form 

P L (n)=A(L)-n- TL , (24) 

limited to the range between n min (L) and n max (L). For simplicity we will assume that the 
lower cutoff is fixed to n m ^{L) = 1. As in the case of the sandpile models, we will assume 
that the upper cutoff is limited by the system size, L, and that n max (L) — > oo when L — > oo. 
The probability distribution Pi{n) should satisfy the normalization condition 

P L (n)dn = 1, (25) 

namely, A(L) = (1 — — 1). In order for A(L) to converge to a finite nonzero 

value as L — > oo, the exponent tl must satisfy > 1 in the infinite system limit. The first 
moment of the distribution, 

(n) L = / nP L (n)dn } (26) 

thus takes the form 

{ ) _ (l~r L )K ax (L) 2 -^-l] 

[ 1L (2-r i )[n max (L) 1 — - 1] ' 1 } 

We observe that for tl > 2 the first moment converges to a finite value in the infinite system 
limit. On the other hand, for 1 < tl < 2, the first moment, (ti)l, diverges for L —>■ oo. We 
thus obtain a connection between the behavior of the first moment of the distribution in the 
infinite system limit and the range of values that the exponent can take. 

For the sandpile models studied here the exponent r L is in the range 1 < t l < 2, and 
indeed, the average avalanche size diverges according to (n) = aL 2 as L — ► oo. The upper 
cutoff n max can be expressed as a function of L and r L : 
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(28) 



" g(2 - r L ) 

_ T L -1 

Using this upper cutoff in the calculation of higher moments we obtain that in the infinite 
system limit they will scale according to Eq. (f2*3*|) with 

2 - r L 

Note that Eq. ()29|) predicts significantly different values of a(q) for the BTW and the Manna 
models, due to the difference in the values of tl for the two models. On the other hand, 
Fig. IHlshows nearly identical values of cr(q), q = 2,3 for the two models. The fact that these 
two moments coincide seems to be due to the deviations from power-law behavior near the 
upper-cutoffs. The effect of these deviations is significant for high moments. 

Recentl y m ultifractal scaling was observed in the avalanche size distribution of the BTW 



model 129 



56 



. This indicates that a finite size scaling analysis of the form of Eq. (|24|) is 



not sufficient for describing the scaling be 
to apply in the case of the Manna model 



ravior of the BTW model, although it was found 



Exponents tl in the range 1 < tl < 2, were observed empirically in the distribution of 
earthquake magnitudes. Many other systems exhibit values of tl in the range 2 < tl < 3. 
In these systems the first moment is kept finite in the infinite system limit, while the second 
moment, that characterizes the fluctuations in the system diverges. Consider, for example, 
a directed graph model describing an internet-like network. Each node in the graph has a 
fixed number, r, of links pointing outwards to other nodes. The graph is constructed such 
that the probability of each node to receive links from newly added nodes is proportional 
to the number of incoming links that it already has. For a network that reached a size of 
L nodes, this process generates a power-law distribution of the number of incoming links 
among the nodes. In the resulting network the total number of outgoing links must be equal 
to the total number of incoming links. Since each node has r outgoing links, the average 
number of incoming links per node must be (n)z, = r, independent of the size L of the 
network. Since the first moment of the distribution is kept finite, while the second moment 
diverges as L —>■ oo, the exponent tl must be in the range 2 < tl < 3 in the infinite system 
limit. The notable feature of the network system is that the average, (n) l, of the power-law 
distribution of the incoming links is forced to remain constant and independent of the system 
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size. Systems that have this feature are common. Other examples include the distribution 
of the number of citations to scientific papers. Each citation is a directed link from a newer 
paper to an older one. While the distribution of the number of outgoing links per paper is 
narrow, the distribution of incoming links is broad, and resembles a power-law distribution. 
Another example is the distributions of income and wealths in western societies, that were 
found to exhibit power-law behavior, at least in the high income sectors, with exponents in 
the range 2 < tl < 3. Here the argument is not as easy to make. However, one may argue 
that the average of these distributions must be connected to the average productivity per 
worker. This productivity remains finite when the size of the economy increases. 

VII. SUMMARY AND CONCLUSIONS 

Abelian sandpile models (both deterministic and stochastic) and random walk models 
have been studied on finite square lattices with open boundaries. The avalanche size distri- 
butions of the sandpile models, as well as the distributions of the lengths of the random-walk 
paths were calculated using various methods. It was shown that, due to the conservation 
laws, the averages (n) of the avalanche size distributions of the deterministic and stochastic 
models are the same, and that they are both equal to the average length of the random-walk 
paths starting from random sites on the same lattice. 
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FIG. 1: The square lattice of size 3 x 3, in which the sites are grouped and indexed according to 
the symmetry of their positions. The indexing starts from the central site and outwards. 

FIG. 2: The average avalanche size for the BTW and Manna models, and the average length of the 
random-walk paths (starting from random sites) on a finite square lattice with open boundaries, 
vs. the lattice size L. For any value of L, it is found that the average avalanche sizes and the 
average path lengths are all equal. 

FIG. 3: Scaling function of the path-length distribution of a random walk on finite lattices of sizes 
L = 32,64,128,256 and 512. Here, n is the walker length, {n)i is its average and P{u)l is the 
probability to obtain a path of length n on a lattice of size L. The five graphs coincide with each 
other. The slope in the linear range is —1/2, as expected for the random walk. 

FIG. 4: The rescaled avalanche size distributions for the BTW model (a) and for the Manna 
model (b) for lattice sizes L = 32, 64, 128, 256 and 512. It is observed that for the BTW model 
the slope exhibits some dependence on L. The slope for best fit obtained by linear regression for 
lattice size 512 corresponds to tl = 1.12 ± 0.02 for the BTW model, and tl = 1-27 ± 0.02 for the 
Manna model. 

FIG. 5: The first three moments of the distribution of the path-lengths of the random walk vs. 
system size. The filled symbols are the results of direct simulations. The empty symbols are the 
results of calculations using the Soler method. The slopes of the best linear fits are 1.98 ±0.04, 
3.97 ±0.06, and 5.96 ±0.08, for q = 1,2 and 3 respectively. 

FIG. 6: The first three moments of the avalanches sizes in the BTW and the Manna models, vs. 
system size. The slopes of the best linear fits for the BTW model are 1.98 ± 0.02, 4.68 ± 0.04, and 
7.52 ± 0.08 for q = 1, 2 and 3 respectively. The slopes of the best linear fits for the Manna model 
are 1.97 ± 0.02, 4.73 ± 0.04, and 7.48 ± 0.08 for q = 1, 2 and 3 respectively. The results for the first 
moment are the same for the two models and coincide with the random walk model. The results 
for higher moments of the two models are not identical, although the differences are small. They 
both are very different from the random walk results. 
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